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Stability analysis and simulations of coupled bulk-surface 

reaction-diffusion systems 
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Abstract. 

In this article we formulate new models for coupled systems of bulk-surface reaction- 
diffusion equations on stationary volumes. The bulk reaction-diffusion equations are cou¬ 
pled to the surface reaction-diffusion equations through linear Robin-type boundary condi¬ 
tions. We then state and prove the necessary conditions for diffusion-driven instability for 
the coupled system. Due to the nature of the coupling between bulk and surface dynamics, 
we are able to decouple the stability analysis of the bulk and surface dynamics. Under a 
suitable choice of model parameter values, the bulk reaction-diffusion system can induce 
patterning on the surface independent of whether the surface reaction-diffusion system 
produces or not, patterning. On the other hand, the surface reaction-diffusion system can 
not generate patterns everywhere in the bulk in the absence of patterning from the bulk 
reaction-diffusion system. For this case, patterns can only be induced in regions close 
to the surface membrane. Various numerical experiments are presented to support our 
theoretical findings. Our most revealing numerical result is that, Robin-type boundary 
conditions seem to introduce a boundary layer coupling the bulk and surface dynamics. 
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1 Introduction 

In many fluid dynamics applications and biological processes, coupled bulk-surface partial 
differential equations naturally arise in (2 D + 3D). In most of these applications and 
processes, morphological instabilities occur through symmetry breaking resulting in the 
formation of heterogeneous distributions of chemical substances m In developmental 
biology, it is essential the emergence and maintenance of polarised states in the form of 
heterogeneous distributions of chemical substances such as proteins and lipids. Examples 
of such processes include (but are not limited to) the formation of buds in yeast cells 
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and cell polarisation in biological cells due to responses to external signals through the 
outer cell membrane [241 [25] . In the context of reaction-diffusion processes, such sym¬ 
metry breaking arises when a uniform steady state, stable in the absence of diffusion, is 
driven unstable when diffusion is present thereby giving rise to the formation of spatially 
inhomogeneous solutions in a process now well-known as the Turing diffusion-driven in¬ 
stability [28] . Classical Turing theory requires that one of the chemical species, typically 
the inhibitor, diffuses much faster than the other, the activator resulting in what is known 
as the long-range inhibition and short-range activation mm- 

Recently, there has been a surge in studies on models that coupled bulk dynamics to 
surface dynamics. For example, Ratz and Roger [25] study symmetry breaking in a bulk- 
surface reaction-diffusion model for signalling networks. In this work, a single diffusion 
partial differential equation (the heat equation) is formulated inside the bulk of a cell, while 
on the cell-surface, a system of two membrane reaction-diffusion equations is formulated. 
The bulk and cell-surface membrane are coupled through Robin-type boundary conditions 
and a flux term for the membrane system [25]. Elliott and Ranner [5] study a finite element 
approach to a sample elliptic problem: a single elliptic partial differential equation is posed 
in the bulk and another is posed on the surface. These are then coupled through Robin- 
type boundary conditions. Novak et al. [[72] present an algorithm for solving a diffusion 
equation on a curved surface coupled to a diffusion model in the volume. Checkkin et al. 
[3] study bulk-mediated diffusion on planar surfaces. Again, diffusion models are posed in 
the bulk and on the surface coupling them through boundary conditions. In the area of 
tissue engineering and regenerative medicine, electrospun membrane are useful in applica¬ 
tions such as filtration systems and sensors for chemical detection. Understanding of the 
fibres’ surface, bulk and architectural properties is crucial to the successful development 
of integrative technology. Nisbet et al. [21] presents a detailed review on surface and bulk 
characterisation of electrospun membranes of porous and fibrous polymer materials. To 
explain the long-range proton translocation along biological mombranes, Medvedev and 
Stuchebrukhov m propose a model that takes into account the coupled bulk-diffusion 
that accompanies the migration of protons on the surface. More recently, Rozada et al., 
[26] presented singular perturbation theory for the stability of localised spot patterns for 
the Brusselator model on the sphere. 

In most of the work above, either elliptic or diffusion models in the bulk have been 
coupled to surface-elliptic or surface-diffusion or surface-reaction-diffusion models posed 
on the surface through Robin-type boundary conditions [31 El [321 ED EMESES]. Here, our 
focus is to couple systems of reaction-diffusion equations posed both in the bulk and on 
the surface, setting a mathematical and computational framework to study more complex 
interactions such as those observed in cell biology, tissue engineering and regenerative 
medicine, developmental biology and biopharmaceuticals [31 El UTl EH E21 EH EH [30]. We 
employ the bulk-surface finite element method as introduced by Elliott and Ranner in [5] to 
numerically solve the coupled system of bulk-surface reaction-diffusion equations. Details 
of the surface-finite element can be found in [4]. The bulk and surface reaction-diffusion 
systems are coupled through Robin-type boundary conditions. The coupled bulk-surface 
finite element algorithm is implemented in deall II [lj. 

The key contributions of our work to the theory of pattern formation are: 

• We derive and prove Turing diffusion-driven instability conditions for a coupled 
system of bulk-surface reaction-diffusion equations. 

• Using a bulk-surface finite element method, we approximate the solution to the 
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model system within the bulk and on the boundary surface of a sphere of radius 
one. 


• Our results show that if the surface-reaction-diffusion system has the long-range 
inhibition, short-range activation form and the bulk-reaction-diffusion system has 
equal diffusion coefficients, then the surface-reaction-diffusion system can induce 
patterns in the bulk close to the surface and no patterns form in the interior, far 
away from the surface. 

• On the other hand, if the bulk-reaction-diffusion system has the long-range inhibi¬ 
tion, short-range activation form and the surface-reaction-diffusion system has equal 
diffusion coefficients, then the bulk-reaction-diffusion system can induce pattern for¬ 
mation on the surface. 


• Furthermore, we prove that if the bulk and surface reaction-diffusion systems have 
equal diffusion coefficients, no patterns form. 

• These theoretical predictions are supported by numerical simulations. 

Hence this article is outlined as follows. In Section [2] we present the coupled bulk- 
surface reaction-diffusion system on stationary volumes with appropriate boundary con¬ 
ditions coupling the bulk and surface partial differential equations. The main results of 


this article are presented in Section 2.2 where we derive Turing diffusion-driven instability 
conditions for the coupled system of bulk-surface reaction-diffusion equations. To vali¬ 
date our theoretical findings, we present bulk-surface finite element numerical solutions in 
Section [3] In Section |4j we conclude and discuss the implications of our findings. 


2 Coupled bulk-surface reaction-diffusion systems on sta¬ 
tionary volumes 

In this section we present a coupled system of bulk-surface reaction-diffusion equations 
(BSRDEs) posed in a three-dimensional volume as well as on the boundary surface enclos¬ 
ing the volume. We impose Robin-type boundary conditions on the bulk reaction-diffusion 
system while no boundary conditions are imposed on the surface reaction-diffusion system 
since the surface is closed. 

2.1 A coupled system of bulk-surface reaction-diffusion equations (BSRDEs) 

Let H be a stationary volume (whose interior is denoted the bulk) enclosed by a compact 
hypersurface F := <9P which is C 2 . Also, let / = [0,T] (T > 0) be some time interval. 
Moreover, let v denote the unit outer normal to T, and let U be any open subset of 
containing T, then for any function u which is differentiable in U, we define the tangential 
gradient on T by, Vr u = Vu — (Vu • u) v, where • denotes the regular dot product and 
V denotes the regular gradient in M Ar+1 . The tangential gradient is the projection of the 
regular gradient onto the tangent plane, thus Vr u ■ i> = 0. The Laplace-Beltrami operator 
on the surface T is then defined to be the tangential divergence of the tangential gradient 
Apu = Vr • Vpu. For a vector function u = (ui, U 2 , ■ ■ ■ , wjv+ 1 ) £ M 7V+1 the tangential 
divergence is defined by 

iV+l 

Vr • u = V • u — ^2 Vi. 

1=1 
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To proceed, we denote by u : 17 x I —> M and two chemical concentrations 

(species) that react and diffuse in 17 and and r : T x I —> 1 and s : T x I —x M be two 
chemical species residing only on the surface T which react and diffuse on the surface. 
In the absence of cross-diffusion and assuming that coupling is only through the reaction 
kinetics, we propose to study the following non-dimensionalised coupled system of BSRDEs 


j ut = V 2 u + 7q/(u,u,), 
Vt = dnV 2 v + ”/ng(u,v), 


in 17 x (0, T], 


n = Vpr + yr (/(r, s) - hi (u, v,r,s)^j, 
s t = dr^ls + 7 r (g(r, s) - h 2 (u, v, r, s)), 


( 2 . 1 ) 


on T x (0, T], 


with coupling boundary conditions 

W ='yrhi(u,v,r,s), 
dn§l =7r h 2 (u,v,r,s), 

In the above, V 2 = represents the Laplacian operator, do, and dr are a 

positive diffusion coefficients in the bulk and on the surface respectively, representing the 
ratio between u and v, and r and s, respectively, jn and yr represent the length scale 
parameters in the bulk and on the surface respectively. In this formulation, we assume that 
/(•, •) and g(-, •) are nonlinear reaction kinetics in the bulk and on the surface. hi(u, v, r, s ) 
and h 2 (u,v,r, s) are reactions representing the coupling of the internal dynamics in the 
bulk 17 to the surface dynamics on the surface T. As a first attempt, we will consider a 
more generalised form of linear coupling of the following nature [12] 


on T x (0, T]. 


( 2 . 2 ) 


hi(u,v,r, s) = air — Piu — Kiv, (2-3) 

h 2 (u,v,r,s) = a 2 s - /3 2 u - k 2 v, (2.4) 


where «i, a 2 , j3i, j3 2 , «i and k 2 are constant non-dimensionalised parameters. Initial 
conditions are given by the positive bounded functions uq(x), vq(x), tq(x) and sq(x). 


2.1.1 Activator-depleted reaction kinetics: An illustrative example 

From now onwards, we restrict our analysis and simulations to the well-known activator- 
depleted substrate reaction model U ES ESI £3 also known as the Brusselator given 
by 

f(u,v)=a — u + u 2 v, and g(u,v) = b — u 2 v, (2-5) 


where a and b are positive parameters. For analytical simplicity, we postulate the model 
system (2.1) in a more compact form given by 

\u t = V 2 u + /i(u, v,r, s), 


= dnV 2 v + f 2 (u, v, r, s ), 

I n = V^r + f 3 (u,v,r,s), 

„ Is* = drV^s + f 4 {u,v,r, s), 


x on 17, t > 0, 


x on T, t > 0, 


( 2 . 6 ) 
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with coupling boundary conditions (2.2)-(2.4). In the above, we have defined appropriately 


fi(u,v,r, s) = 70 (a - u + u 2 v), 
f 2 (u,v,r,s) = ^(b-u 2 v), 

f 3 (u,v,r, s ) = 7r(o - r + r 2 s - a 3 r + fiiu + k 3 v), 
/4(u, v, r, s ) = 7 r (b — r 2 s — a 2 s + @2u + n 2 v) . 


(2.7) 

( 2 . 8 ) 
(2.9) 

( 2 . 10 ) 


2.2 Linear stability analysis of the coupled system of BSRDEs 

Definition 2.1 (Uniform steady state). A point (it*, v*,r*, s*) is a uniform steady state of 
the coupled system of BSRDEs (2.6) with reaction kinetics (|2.5l) if it solves the nonlinear 


algebraic system given by fi(u*,v*,r *, s*) = 0, for all i = 1, 2, 3, 4, and satisfies the 


boundary conditions given by (2.2)-(2.4). 


Proposition 2.1 (Existence and uniqueness of the uniform steady state). The coupled 
system of BSRDEs ( |2.6| ) with boundary conditions (2.2) admits a unique steady state given 
by 

(»*,»*,r*, 8 *) = (a + b.Df^.a + b.D^) , (2.11) 

provided the following compatibility condition on the coefficients of the coupling is satisfied 

(fii — au)(« 2 — ot 2 ) — fi 2 = 0. (2.12) 

Proof. The proof follows immediately from the definition of the uniform steady state 


satisfying reaction kinetics (2.7)-(2.10). It must be noted that in deriving this unique 

ics 

□ 


uniform steady state the compatibility condition ( 2 . 12 ) coupling bulk and surface dynamics 
must be satisfied. 


Remark 2.1. The constraint condition (2.12) on the parameter values on, ffi, and k *, i = 1, 2 
is a general case of the specific parameter values given in |T2] where the following parameter 
values where selected qi = fi\ = -jf, a 2 = n 2 = 5, k\ = 0, and fi 2 = 0. 


Remark 2.2. Note that there exists an infinite number of solutions to problem (2.12). 


2.2.1 Linear stability analysis in the absence of diffusion 


Next, we study Turing diffusion-driven instability for the coupled system of BSRDEs (2.1)- 


(2.4) with reaction kinetics (2.5). To proceed, we first consider the linear stability of the 
spatially uniform steady state. For convenience’s sake, let us denote by w = (u, v, r, s') , 
the vector of the species it, v, r and s. Furthermore, defining the vector £ such that 
|£i| << 1 for all i = 1, 2, 3 and 4, it follows that writing w = w* +£, the linearized system 
of coupled BSRDEs can be posed as 


w t = = J f £> 


(2.13) 
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where J F represents the Jacobian matrix representing the first linear terms of the lin¬ 
earization process. Its entries are defined by 


J ]? = 


/dfi 

dii 

ah 

& 

du 


SA 

ii 

dv 

dh 

$L 


dh 

dr 

dh 

dr 

dfs 

dr 

dh 

dr 

dfi\ 
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ds 


(flu fly 

f^U f^v 

0 

0 

0 0 



oh 

ds 


fSu f3v 

hr 

h s 



dfi 

ds / 


f 4 u 

hr 

hs) 






( fu 

fv 

0 


0 ^ 




9u 

9v 

0 


0 




hi u 

~h\ v 

fr- 

h\ r fs hi s 




\ ^2 u 

~h j 2y 

9r — 

h2 r 9s ~ h2 s J 


(2.14) 


where by definition f\ u := represents a partial derivative of with respect to u. 

We are looking for solutions to the system of linear ordinary differential equations (2.13) 
which are of the form £ oc e At . Substituting into (2.13), results in the following classical 
eigenvalue problem 


XI-J, 


= 0 , 


(2.15) 


where I is the identity matrix. Making appropriate substitutions and carrying out stan¬ 
dard calculations we obtain the following dispersion relation for A 


XI- J 


^ — hu hv 0 0 

h u ^ ~ h v 0 0 

hu hv ^ - hr hs 

f4 u .f4v f4 r ^ ./J 

Pi( A) A‘ OqA^ -(- CL 2A^ T U3A (T4 = 0 , 


= 0, 


where 


®i — — (flu + h v + hr + /4s) > 

«2 = ( fluhv ~ flvhu) + (hrhs - hshr ) + (hu + hv)(hr + /4s)) 


«3 = 


(fluhv - flvhu)(hr + hs) + (hrhs - hshr)(hu + hi 
a 4 = (fluhv - flvhu)(hrhs ~ hshr)- 

For convenience’s sake, let us denote by 


( J F)n : = 


flu hv 

h u h v 


and (J p) r : = 


hr hs 
fir fis 


(2.16) 

(2.17) 

(2.18) 
(2.19) 


( 2 . 20 ) 


the submatrices of Jp corresponding to the bulk reaction kinetics and the surface reaction 
kinetics respectively. We can now define 


Tr (J F ) ■■= flu + hv +hr + hs> Tl '( J F)n : =hu + h v , ^ {J f)t := h r + hs, 
Det (J F ) n := h u h v ~ h v hu , and Det ( Jf)t := hrhs ~ hshr- 

Theorem 2.1 (Necessary and sufficient conditions for Re(A) < 0). The necessary and 
sufficient conditions such that the zeros of the polynomial Pi(X) have Re( A) < 0 are given 
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by the following conditions 


Tr {J F ) < 0, 

Det (J j p)q + Det (J F ) r + Tr {J p)n Tr (J F ) r > 0, 

Det (J F )q Tr (Jp) r + Det (J p) r Tr (J p)n < 0, 

Det ( J p)n Det (Jp) r > 0, 

Tr ( j f)t Tr i J f) ~ 2Det( K J F)n] Tr ( J F)n + [ Tr ( J F)n Tr ( J f) 


( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 

(2.24) 


— 2Det (J F ) r 


Tr {Jp)r > °> 


{Det (J f)m + Det ( J f)v ) {Det {J f)q Dr {Jp)r 

+Det {J F )r Tr ( J F)n) Tr ( J F )] Tr ( J f)o Tr ( J F)r > °- 


(2.25) 


(2.26) 


Proof. The proof enforces that 734 (A) is a Hurwitz polynomial and therefore satisfies the 
Routh-Hurwitz criterion in order for Re(A) < 0. The first condition to be satisfied is that 
«4 7 ^ 0 otherwise A = 0 is a trivial root, thereby reducing the 4-th order polynomial to a 
cubic polynomial. The first four conditions are a result of requiring that each coefficient a* 
with i = 1, 2, 3 and 4 of the polynomial 734 (A) are all positive. The rest of the conditions 
are derived as shown below. 

We require that the determinant of the matrix 


01 03 . n 

— aia 2 — 03 > (J. 

1 a 2 

Substituting a±, a 2 and 03 appropriately we obtain 

Det (■ J F)n + Det (■ J F)r + Tr (■ J F)n Tr (• J f)t } [ “ Tr (■ J f) 
+ [Det {J F ) n Tr (J F ) r + Det (' J f)t ^ (■ J f)q 


> 0. 


Exploiting the fact that 


Tr(J F ) =Tr(J F ) n + Tr(J F ) r , 


(2.27) 


it then follows that 


aia 2 — 03 

= Tr ( J f)q T i ' ( j f) r Tr ( J f) “ Det (■ J F)n Tr ( J f)q + Det (■ J f)t Tr ( J f)i 

if and only if 
1 IT 


> 0 


2 L 


Tr (Jp)pTr {j F ) ~ 2 Det (Jp) n J Tr («Jp)n 
1 


+ 2 


Tr ( Tr (■ J f) ~ 2Det (■ J f)t Tr (■ J f)v > °- 


Multiplying throughout by 2 we obtain condition ( |2.25 ) in Theorem |2.1 
The last condition results from imposing the condition that 


a \ 

03 

0 

1 

02 

04 

0 

Oi 

03 


03(0.1(12 — 03) — 0^04 > 0. 
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It can be shown that 


03(0102 — 03) 

= - Det (J F ) n Tr (Jp) n Tr 2 (J F ) r + Det (J F ) r Tr 2 ( J F ) n Tr (J F ) r ] Tr (Jp) 

+ Det 2 (J F ) 0 Tr (J F ) n Tr (J F ) r + Det (J F ) n Det (J F ) r Tr 2 (J F ) n 
+ Det (J F ) n Det (J F ) r Tr 2 (J F ) n + Det 2 (J F ) r Tr (J F ) n Tr (J F ) r . (2.28) 

On the other hand, 

a 2 04 = Tr 2 (Jp) Det (J F )q Det (J F ) r 
= (Tr(J F ) n + Tr(J F ) r ) Det (Jp) Q Det (Jp) r 

= Det (J F ) n Det (J F ) r Tr 2 (J F ) n + 2Det (J F ) n Det (Jp) r Tr (J F ) n TV (J F ) r 

+ Det (J F ) n Det (J F ) r Tr 2 (J F ) r • (2.29) 

Hence combining (2.28) and ( |2.29[ ) and simplifying conveniently we have 

03(0102 - o 3 ) - a 2 a 4 = (Det (J F ) n + Det (J F ) r ) 2 - (Det (J F )^Tr (J F ) r 
+Det (Jp) r Tr(J F ) n )Tr(J F )] x Tr (J F ) n Tr (J F ) r >0, 
resulting in condition ( 2.26[ ). □ 

Remark 2.3. The characteristic polynomial 

p 4 (A) = A 4 + aiA^ + <22 A 2 + (X3 A + a 4 

can also be written more compactly in the form of 

Pa{ A) = ^A“ + A(/i u + f2 v ) + fl u f2v — flvhu) + ^(/3 r + f^s) + hrf^s ~ 

thereby coupling the bulk and surface dispersion relations in the absence of spatial varia¬ 
tions. 


2.2.2 Linear stability analysis in the presence of diffusion 

Next we introduce spatial variations and study under what conditions the uniform steady 
state is linearly unstable. We linearize around the uniform steady state by taking small 
spatially varying perturbations of the form 


w(x 

,t) = w* + e£(x,t), with e << 1. 

(2.30) 

Substituting (2.30) into the coupled system of BSRDEs (2.1)-(2.4) with reaction kinetics 
(2.5) we obtain a linearized system of partial differential equations 

Ci t = v2 £i + 7o 

/u£l + fv £ 2 ^, 

(2.31) 

£2 1 = dnV 2 ^ + 7n {duii + 9vi 2 ), 

(2.32) 

£ 3 1 = V 2 £3 + yr 

/r£ 3 + /s£4 - ^lu£l - /llu£2 “ ^lr£3 ~ , 

(2.33) 

£44 = drVp£ 4 + 7r(gr£3 + 9sU — /i2 u £i — h 2v C 2 - ^£3 ~ ^£ 4 ) ■ 

(2.34) 
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with linearised boundary conditions 

= 7r (hiu€i + hi v %2 + hi r £ 3 + /ii s £ 4 ) , (2.35) 

= 7r (h 2 u £,i + h r 2 v & + for&s + h 2 s t,±J. (2.36) 

In the above, we have used the original reaction kinetics for the purpose of clarity. 

In order to proceed, we restrict our analysis to circular and spherical domains where 

we can transform the cartesian coordinates into polar coordinates and be able to exploit 
the method of separation of variables. Without loss of generality, we write the following 
eigenvalue problem in the bulk 


v V*,, m (r) = -klm^k liTn {r), 0 < r < 1 , 


(2.37) 


where each ip k satisfies the boundary conditions (2.35) and (2.36). On the surface the 
eigenvalue problem is posed as 


Vr^(y) = -l(l + i)<Ky), y e r. 


(2.38) 


Remark 2.4. For the case of circular and spherical domains, if r = 1, then kf m = l(l + 1). 

Taking x E B, y E T, then writing in polar coordinates x = ry , r E (0,1) we can define, 
for all l E Nq, m£Z, \m\ < l, the following power series solutions |24l [25] 


£i(ry,t) = ^ui t m(t)i>k l , rn (r)<pi,m(y), &(ry,t) =^^, m (l)4,, m (’')^,m(y), (2.39) 

£3(y?^) ^ ^ and £4(2/?^) E Sl,m (t)(pl,m(y)- (2.40) 


On the surface, substituting the power series solutions (2.40) into (2.33) and (2.34) we 
have 


dri 


dt = -1(1 + 1 )ri t m + 7r [frn,m + fsSl , m ) 

7r(hiuUi i m'tl>k l m (l) h / ±v'Vl,m'lfiki m (1) d - ^1 r^l,m d - ^1 s$l,r 


ds t . 


dt = -drl(l + l)si t m+'Yr(9rn,m +gsSl,m^ 

~ 7r( h2 u ui,m'ipki :rn (l) + h2 v v^ m ip klm (l) + /i 2r n jm + h 2s si, r 


(2.41) 


(2.42) 


Similarly, substituting the power series solutions (2.39) into the bulk equations (2.31) and 
(2.32) we obtain the following system of ordinary differential equations 

duijm 


dt 

dvi,n 

dt 


— ^ 1 , 171 ^ 1 ,Tn T 7fi (^fuUl^m T • 

diik Lrn V Lm + 7f2 (^fJuMl.m + • 


(2.43) 

(2.44) 

(2.45) 

dnvi, m ip' kl ( 1) = 7 r(/i 2 u ^,mV , fc i , m (l) + h 2v vi jrn 'ip klm (l) + h 2r .n,m + h 2s si )m V (2.46) 


Equations (2.43) and (2.44) are supplemented with boundary conditions 

m (1) = 7r ( h\ U ui tm il>ki m (1) T hi v vi : m' l Pki m (1) T hi r ri, m + h\ s si jT 
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where M := 


dr 


r =1 


Writing 


o \ 


( u l,mi v l,mi r l,mi s l,mi) ~ ( u l,m’ V l,m’ r l,m’ S l,m) e 


and substituting into the system of ordinary differential equations (2.41)-(2.44), we obtain 
the following eigenvalue problem 


(\ l , m I + M)tf tm = 0 


(2.47) 


where 


M = 


( k lm - 7 nfu 

~'yo.9u 

V 0 


and 


-7 nfv 0 0 \ 

( h k l m - ingv 0 0 

0 1(1 + 1) — 7r/r —7r/s 

~lv9r d r Z(Z + 1) - 7r 9s) 

fiO = /„.o o „o „0 \ T 

S <l,m \ u l,m' u l,mi l,mi 0 l,m) 


Note that the boundary conditions ( |2.45 ) and (2.46) have been applied appropriately to 
the surface linearised reaction-diffusion equations. Since 

('C C 4’ + (°> °- °> °) T > 

it follows that the coefficient matrix must be singular, hence we require that 


^Lm,I + M 


= 0. 


Straight forward calculations show that the eigenvalue A i m solves the following dispersion 
relation written in compact form as 

+ Tr (M) a A i tTn + Det (M) a ^ ^A f m + Tr (M) r A+ Det (M) r j = 0, (2.48) 

where we have defined conveniently 

Tr(M)a := (dn + l)& 2 m - 7 n(/« + S«), 

Tr(M) r := (d r + 1)Z(* + 1) - 7r (fr + 9s), 

Det(M)a := - 70 (dn/ u + g v ) fc 2 m + 7 - fv9u), 

Det (M) r := d r ^ 2 (^ + l) 2 - 7r (dr/r + &) ^ + 1) + 7r (fr9s ~ fs9r)- 

The above holds true if and only if either 

rf,m + Tr (M) n \m + Det (M) n = 0, (2.49) 

or 

A 2 m + Tr (l\T) r A i }m + Det (AT) r = 0. (2.50) 

In the presence of diffusion, we require the emergence of spatial growth. In order for the 
uniform steady state w* to be unstable we require that either 

1. Re(Xi jrn (kf m )) > 0 for some kf m > 0, 
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or 


2. Re(A+ 1))) > 0 for some l(l + 1) > 0, 


or 

3. both. 


Solving (2.49) (and similarly (2.50)) we obtain the eigenvalues 


2Re(A i,m(klJ) = -Tr (M) n ± ^Th 2 (M) n - 4Det (M) n . (2.51) 

It follows then that Re(A i, m (kf m )) > 0 for some kf m > 0 if and only if the following 
conditions hold: 


or 


Tr (M) n < 0 (dn + l)kf m - 70 (f u + g v ) < 0, and 
„Det (M) n > 0 <=^> dnkf m - 70 {d n f u + g v ) k 2 m + 7 ^(f u g v - f v 9u) > 0, 

Tr (M) n > 0 (do + l)fe 2 m - 7 n(/ u + £„) > 0, and 

„°et (M) n < 0 4=> dnkf m - 70 ( dnfu + £„) + 7 - fvg u ) < 0. 


(2.52) 


(2.53) 


Similarly, on the surface, Re(A i t m(l(l + 1))) > 0 for some /(I + 1) > 0 if and only the 
following conditions hold: 


Tr(M) r <0 4=> (d r + 1)1(1 + 1) - 7r(/r + g s ) < 0, and 

< 

,Det (M) r > 0 d r l 2 (l + l) 2 - 7r (dr/r + 9s) l (l + 1) + 7r (fr9s ~ fs9r) > 0, 

(2.54) 
or 

Tr(M) r >0 4=> (d r + 1)1(1 + 1) — 7r(/r + gs) > 0, and 

< 

,Det (M) r < 0 <==> d r l 2 (l + l) 2 - 7r (dr/r + g s ) l (l + 1) + 7r(/r5s - fs9r) < 0. 

(2.55) 

We are in a position to state the following theorem: 

Theorem 2.2. Assuming that 

Tr i J F)n = /« + 9v < 0 and Det (J F ) n = f u g v - f v g u > 0, (2.56) 

then the necessary conditions for Re(\i tTn (kf m )) > 0 for some k tm > 0 are given by 

dnfu + 9v > 0, and (dnfu + 9vf - 4 dn ( f u g v - fv9u ) > 0. (2.57) 

Similarly, assuming that 

Tr(j f) t = fr + g s and Det(J F ) r = f r g s - f s g r > 0, (2.58) 

then the necessary conditions for Re(X^ m (l(l + 1))) > 0 for some l (l + 1) > 0 are given by 
dr/r + 9s > 0, and (d r /r + gs) 2 ~ 4 d r (f r g s ~ fs9r) > 0. (2.59) 
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Proof. The proof is a direct consequence of conditions (|2.52 )- (2.55). Assuming conditions 
(2.56) and (2.58) hold, then one of the conditions in (2.52) and ( |2.54 ) is violated, which 
implies that Re(Az,m(^ 2 m )) < 0 for all kf m > 0 and similarly Re(A i : m(l(l + 1))) < 0 for all 
1(1 + 1) >0. This entails that the system can no longer exhibit spatially inhomogeneous 
solutions. 

The only two conditions left to hold true are (2.53) and (2.55). This case corresponds 
to the classical standard two-component reaction-diffusion system which requires that (for 
details see for example USD 


dnfu +g v > 0, and (dnf u + 9vY ~ 4 do. ( f u g v - f v g u ) > 0, 

and similarly 

drfr + 9s> 0, and (d r f r + g s Y - 4 d r ( f r g s ~ fs9r) > 0. 
This completes the proof. 


(2.60) 

(2.61) 

□ 


Remark 2.5. Assuming conditions (2.56) and (2.58) both hold, then conditions (2.57) and 
(2.59) imply that dn Y 1 and dr Y 1- 


Remark 2.6. If condition (2.56) or (2.58) holds only, then either Y 1 or dr Y 1 but not 
necessarily both. 


Remark 2.7. If conditions (2.56) and (2.58) are both violated, then diffusion-driven insta¬ 
bility can not occur. 


Remark 2.8. Similarly to classical reaction-diffusion systems, conditions (2.57) and (2.59) 
imply the existence of critical diffusion coefficients in the bulk and on the surface whereby 
the uniform states lose stability. In order for diffusion-driven instability to occur, the bulk 
and surface diffusion coefficients must be greater than the values of the critical diffusion 
coefficients. 

Next we investigate under what assumptions on the reaction-kinetics do conditions 


(2.52) and (2.54) hold true. 

• First let us consider the case when 

Tr ( J f)q = fu + 9 v> 0 and Det (Jp) a = fu 9 v ~ fv 9 u > 0, 

and 

Tr (Jp) r = fr + 9 s > 0 and Det (Jp) r = fr 9 s ~ fs 9 r > 0. 


Then Tr (J p) = Tr (Jf)q + Tr («7j i ) r > 0 which violates condition (2.21). 
Similarly the case when 

Tr (■ J F)n = fu + 9v > 0 and Det (Jp) n = f u 9v ~ fv9u < 0, 

and 

Tr (Jp) r = fr + g s > 0 and Det (j^) r = f r g s ~ f s g r < 0 
violates condition (|2.21[). 
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• Let us consider the case when 


Tr ( J f)q = fu + 9v < 0 and Det (Jp) n = f u g v ~ fvQu < 0, 


and 

Tr ( j f)v = f r + 9s < 0 and Det (Jp) r = f r g s - f s g r < 0. 

Then it follows that condition ( j2.25| ) given by 

Tr ( J F)r Tr ( J f) “ 2Det ( j f)o Ti '( J f)q+ ^ ( J J p)n Tr ( J f) 

— 2Det (J p) r Tr (j^) r < 0, 

is violated. 

Next we consider the case when 

Tr ( J F)n = fu + 9v < 0 and Det (Jp) n = f^9v ~ fv9u < 0, 

and 

Tr ( j f)v = f' r + 9s > 0 and Det (Jp) r = f r g s - f s g r < 0. 

It follows then that none of the conditions ( |2.21[ )-( |2.26 ) are violated. However, 
condition (2.52) does not hold. 

Similarly the case when 

Tr ( J F)n = fu + 9v> 0 and Det (Jp) n = f u g v ~ fv9u < 0, 

and 

Tr (Jf) v = fr + 9s < 0 and Det (Jp) r = f r g s - f s g r < 0. 

This implies that that none of the conditions ( |2.21[ )- (2.26) are violated, while con¬ 
dition (2.54) fails not hold. 

Finally, the cases when 


Tr (Jp) n = f u + 9v > 0 and Det (J F ) Q = f u g v - f v g u > 0, 
Tr (Jp) T = fr + 9s < 0 and Det (Jp) r = f r g s ~ f s g r > 0, 


(2.62) 


and 


j Tr (Jp) n = fu + 9v < 0 and Det (Jp) n = fu9v ~ fv9u > 0, 
[Tr (Jp) r = f r + g s > 0 and Det (Jp) T = fr9s - fs9r > 0, 

result in Remark 12.61 above. 


(2.63) 


The above cases clearly eliminate conditions (2.52) and (2.54) as necessary for uniform 
steady state to be driven unstable in the presence of diffusion. We are now in a position 
to state our main result. 
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Theorem 2.3 (Necessary conditions for diffusion-driven instability). The necessary con¬ 


ditions for diffusion-driven instability for the coupled system of BSRDEs (2.1) - (2.4) are 
given by 

Tr {J F ) < 0 , (2.64) 

Det (J F ) n + Det(J F ) r + Tr(j F ) n Tr(j F ) r > 0, (2.65) 

Det(J F ) n Tr(J F ) r + Det(j F ) r Tr(j F ) n < 0, (2.66) 

Det ( J F ) n Det (J F ) r > 0, (2.67) 

Tr ( J F ) r Tr { J F)- 2 Det ( J F)n Tr ( J F)n+ Tr ( J f) n Tr ( J f) 

-2Det(J F ) r ]Tr(J F ) r >0, (2.68) 

(Det(J F ) Q + Det(J F ) r f - (. Det(J F ) Q Tr(J F ) r 

+Det(J F ) r Tr(J F ) n )Tr(J F )] Tr ( J F ) n Tr (Jp) r > 0, (2.69) 


and 

dnfu + 9v > 0, 

and 

{dnfu + 9vf ~ 4do-Def (J F ) n > 0, 

(2.70) 

or/and 

dvfr + g s > 0, 

and 

(dr.fr + 9s) 2 ~ 4drDet (J F ) r > 0. 

(2.71) 

2.2.3 

Theoretical predictions 




From the analytical results we state the following theoretical predictions to be validated 
through the use of numerical simulations. 

1. The bulk and surface diffusion coefficients do and dr must be greater than one in 
order for diffusion-driven instability to occur. Taking do = dr = 1 results in a 


contradiction between conditions (2.64), (2.70) and (2.71). As a result, the BSRDEs 


does not give rise to the formation of spatial structure. For this case, the uniform 


steady state is the only stable solution for the coupled system of BSRDEs (2.1) 


(2.4). 


2. The above imply that taking do > 1 and dr = 1, the bulk-reaction-diffusion system 
has the potential to induce patterning in the bulk for appropriate diffusion-driven 
instability parameter values while the surface-reaction-diffusion system is not able 


to generate patterns. Here all the conditions (2.64)-(2.70) hold except (2.71). 


3. Alternatively taking do = 1 and dr > 1, the bulk-reaction-diffusion system fails 
to induce patterning in the bulk while the surface-reaction-diffusion system has the 


potential to induce patterning on the surface. Similarly, all the conditions (2.64) 


(2.71) hold except (2.70). 


4. On the other hand, taking do > 1 and dr > 1 appropriately, then the coupled 
system of BSRDEs exhibits patterning both in the bulk and on the surface. All the 


conditions (2.64)-(2.71) hold both in the bulk and on the surface. 
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3 Numerical simulations of the coupled system of bulk-surface 
reaction-diffusion equations (BSRDEs) 


In this section we present bulk-surface finite element numerical solutions corresponding 
to the coupled system of bulk-surface reaction diffusion equations (BSRDEs) (2.1)-(2.5). 
Here we omit the details of the bulk-surface finite element method as these are given else¬ 
where (see [15] for details). Our method is inspired by the work of Elliott and Ranner [5]. 
We use the bulk-surface finite element method to discretise in space with piecewise bilin¬ 
ear elements and an implicit second order fractional-step 0-scheme to discretise in time 
using the Newton’s method for the lineraisation ii. For details on the convergence 
and stability of the fully implicit time-stepping fractional-step 0-scheme, the reader is re¬ 
ferred to Madzvamuse et al. mm- In all our numerical experiments, we fix the kinetic 
model parameter values a = 0.1 and b = 0.9 since these satisfy the Turing diffusion-driven 
instability conditions ( |2.64 )-(2.71). For these model parameter values the equilibrium val¬ 
ues are (u*,v*,r*, s*) = (1, 0.9,1, 0.9). Initial conditions are prescribed as small random 
perturbations around the equilibrium values. For illustrative purposes let us take the pa¬ 
rameter values describing the boundary conditions as shown in Table [l] these are selected 
such that they satisfy the compatibility condition (2.12). 


a 

b 

7n 

7r 

a\ 

«2 

Pi 

p2 

Kl 

K2 

0.1 

0.9 

500 

500 

5 

12 

5 

5 

12 

0 

0 

5 


Table 1: Model parameter values for the coupled system of BSRDEs (2.1) - (2.4). 


3.0.4 A note on the bulk-surface triangulation 

We briefly outline how the bulk-surface triangulation is generated. For further specific 
details please see reference |15j . Let Qh be a triangulation of the bulk geometry D with 
vertices Xi, i = 1, ...,Nh, where Nh is the number of vertices in the triangulation. From 
fl/j we then construct T^ to be the triangulation of the surface geometry T by defining 
T/j = Qh\dn h , be. the vertices of T/,, are the same as those lying on the surface of Qh- In 
particular, then, we have dflh = T^. An example mesh is shown in Figure [l] The bulk 
triangulation induces the surface triangulation as illustrated. 


3.1 Numerical experiments 


In this section we will only present four cases to validate our theoretical predictions out¬ 
lined in Section 2.2.3 In most of our simulations parameter values are fixed as shown 


in Table [l] except for dn and dr whose values are varied to demonstrate the patterning 
mechanism of the coupled system of BSRDEs. We only present patterns corresponding 
to the chemical species u and r in the bulk and on the surface respectively. Those corre¬ 
sponding to v and s are 180 degrees out of phase to those of u and r and are therefore 
omitted. It must be noted however that this is not always the case in general, Robin-type 
boundary conditions may alter the structure of the solution profiles depending on the 
model parameter values and the coupling compatibility boundary parameters. 


15 




























Figure 1: Example meshes for the bulk (top) and surface system (bottom). Part of the 
domain has been cut away and shown on the right to reveal some internal mesh structure. 



Figure 2: Numerical solutions corresponding to the coupled system of BSRDEs (2.1)-(2.5) 
with do = 1 in the bulk and dp = 1 on the surface. The uniform steady state solutions are 
converged to and no patterns form. Columns 1 and 2: solutions in the bulk representing 
u. Columns 3 and 4: solutions on the surface representing r. Second and fourth columns 
represent cross sections of the bulk and the surface respectively (Colour figure online). 


3.1.1 Simulations of the coupled system of BSRDEs with (do,dp) = (1,1) 

The bulk-surface finite element numerical simulations of the coupled system of BSRDEs 
with with do = 1 in the bulk, dp = 1 on the surface are shown in Figure [2j We observe 
that no patterns form in complete agreement with theoretical predictions. Similarly to 
classical reaction-diffusion systems, diffusion coefficients must be greater than one. In 
particular, the diffusion coefficients must be greater than their corresponding respective 
critical diffusion coefficients in the bulk and on the surface. An example is shown next. 

3.1.2 Simulations of the coupled system of BSRDEs with (do, dp) = (1,20) 

For illustrative purposes, let us take do = 1 in the bulk, dp = 20 > df" = 8.5 on the 
surface. Figure [3] illustrate pattern formation on the surface as well as within a small 
region in the vicinity of the surface membrane. Spots are observed to form on the surface, 
while in the bulk, small balls form inside. Far away from the surface, no patterns form 
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Figure 3: Numerical solutions corresponding to the coupled system of BSRDEs (2.1)-(2.5) 
with dpi = 1 in the bulk and dp = 20 on the surface. Columns 1 and 2: solutions in the 
bulk representing u. Columns 3 and 4: solutions on the surface representing r. Second and 
fourth columns represent cross sections of the bulk and the surface respectively (Colour 
figure online). Spot patterns form on the surface while small balls form in the vicinity of 
the surface inside the bulk. 


since the necessary conditions for diffusion-driven instability are not fulfilled in the bulk. 
These results confirm our theoretical predictions. We note that this particular example 
describes realistically pattern formation in biological systems. We expect skin patterning 
to manifest in the epidermis layer as well as on the surface. 

3.1.3 Simulations of the coupled system of BSRDEs with (dn, dp) = (20,1) 



Figure 4: Numerical solutions corresponding to the coupled system of BSRDEs (2.1)-(2.5) 
with dpi = 20 in the bulk and dp = 1 on the surface. Columns 1 and 2: solutions in the 
bulk representing u. Columns 3 and 4: solutions on the surface representing r. Second and 
fourth columns represent cross sections of the bulk and the surface respectively (Colour 
figure online). Spectacular patterning occurs in the bulk exhibiting spots, stripes and 
circular patterns. The surface dynamics produce uniform patterning. 


To generate patterns in the bulk we take dpi = 20 > df rit =8.5 and dr = 1 on the 
surface. Figure [4] exhibits stripe, circular and spot patterns in the bulk as illustrated 
by the cross-sections. On the surface, uniform patterns occur consistent with theoretical 
predictions. Although the patterns for the u species (columns one and two) appear uniform 
on the surface this is simply due to the colour scale, with the amplitude of the patterns in 
the bulk larger than those on the surface. This difference in the amplitude of the pattern 
of the bulk solution in the bulk and on the surface is due to the Robin type boundary 
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conditions. Unlike zero-flux (also known as homogeneous Neumann) boundary conditions 
for standard reaction-diffusion systems which imply that no species enter or leave the 
domain, here, there is deposition or removal of chemical species through the flux on the 
surface, resulting in differences in amplitude between the bulk and surface solutions. 

3.1.4 Simulations of the coupled system of BSRDEs with (dn,dr) = (20,20) 



Figure 5: Numerical solutions corresponding to the coupled system of BSRDEs (2.1)-(2.5) 
with dn = 20 in the bulk and dr = 20 on the surface. Columns 1 and 2: solutions in the 
bulk representing u. Columns 3 and 4: solutions on the surface representing r. Second and 
fourth columns represent cross sections of the bulk and the surface respectively (Colour 
figure online). We observe spot pattern formation both in the bulk and on the surface. 


In this example, we illustrate how both bulk and surface dynamics induce patterning 
by taking dn = 20 in the bulk, dr = 20 on the surface. Figure [5] shows pattern formation 
in the bulk and on the surface. In the bulk we observe the formation of balls (which can 
be seen as spots through cross-sections) and these translate to spots on the surface. The 
surface dynamics themselves induce spot pattern formation. 


4 Conclusion, discussion and future research challenges 

We have presented a coupled system of bulk-surface reaction-diffusion equations whereby 
the bulk and surface reaction-diffusion systems are coupled through Robin-type boundary 
conditions. Nonlinear reaction-kinetics are considered in the bulk and on the surface and 
for illustrative purposes, the activator-depleted model was selected since it has a unique 
positive steady state. By using linear stability theory close to the bifurcation point, we 
state and prove a generalisation of the necessary conditions for Turing diffusion-driven 
instability for the coupled system of BSRDEs. Our most revealing result is that the bulk 
reaction-diffusion system has the capability of inducing patterning (under appropriate 
model and compatibility parameter values) for the surface reaction-diffusion model. On the 
other hand, the surface reaction-diffusion is not capable of inducing patterning everywhere 
in the bulk; patterns can only be induced in regions close to the surface membrane. For 
skin pattern formation, this example is consistent with the observation that patterns will 
form on the surface as well as within the epidermis layer close to the surface. We do not 
expect patterning to form everywhere in the body of the animals. 

Our studies reveal the following observations and research questions still to addressed: 

• Our numerical experiments reveal that the Robin-type boundary conditions seem 


18 





to introduce a boundary layer coupling the bulk and surface dynamics. However, 
these boundary conditions do not appear explicitly in the conditions for diffusion- 
driven instability and this makes it difficult to theoretically analyse their role and 
implications to pattern formation. Further studies are required to understand the 
role of these boundary conditions as well as the size of the boundary layer. 


The compatibility condition (2.12) implies that the uniform steady state in the bulk 
is identical to the uniform state on the surface. We are currently studying the 
implications of relaxing the compatibility condition. 


• Finally, in this manuscript, we have not carried out detailed parameter search and 
estimation to deduce the necessary and sufficient conditions for pattern generation 
as well isolating excitable wavenumbers in the bulk and on the surface. Such studies 
might reveal more interesting properties of the coupled bulk-surface model and this 
forms part of our current studies. 

We have presented a framework that couples bulk dynamics (3D) to surface dynamics 
(2D) with the potential of numerous applications in cell motility, developmental biol¬ 
ogy, tissue engineering and regenerative medicine and biopharmaceutical where reaction- 
diffusion type models are routinely used OEunusiEaiMusaEO]. 

We have restricted our studies to stationary volumes. In most cases, biological surfaces 
are known to evolve continuously with time. This introduces extra complexities to the 
modelling, analysis and simulation of coupled systems of bulk-surface reaction-diffusion 
equations. In order to consider evolving bulk-surface partial differential equations, evolu¬ 
tion laws (geometric) should be formulated describing how the bulk and surface evolve. 
Here, it is important to consider specific experimental settings that allow for detailed 
knowledge of properties (biomechanical) and processes (biochemical) involved in the bulk- 
surface evolution. Such a framework will allows us to study 3D cell migration in the area 
of cell motility B 0 El ED]. In future studies, we propose to develop a 3D integrative 
model that couples bulk and surface dynamics during growth development or movement. 
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